program abadiels, eclass
	version 13
	local version 
	/* 
	This program implements Abadie-weighted least squares.
	It uses Mata because Stata won't accomodate negative weights
	The syntax is abadiels y d z x1 x2 ...
	*/
	syntax varlist(min=2 numeric) [if] [in] 
	marksample touse
	
	* initialize temporary names, variables, and files:
	tempname beta varcov ses results
	tempvar kappa
	tempfile 
	tokenize `varlist'
	local y "`1'"
	macro shift
	local d "`1'"
	macro shift
	local z "`1'"
	macro shift
	local x "`*'"
	preserve
	qui keep if `touse'
	
	* generate kappa
	qui sum `z'
	local tau=r(mean)
	gen `kappa'=1-`d'*(1-`z')/(1-`tau')-(1-`d')*`z'/`tau'
	mata output=wls("`y' `kappa' `d' `x'")
	mata beta=*output[1]
	mata varcov=*output[2]
	mata st_matrix("`beta'",beta')
	mata st_matrix("`varcov'",varcov)
	mat colnames `beta'= `d' `x' constant
	mat rownames `beta'= coeff
	mat rownames `varcov'= `d' `x' constant
	mat colnames `varcov'= `d' `x' constant
	mat diag=vecdiag(`varcov')
	mat `ses'=vecdiag(`varcov')
	mata st_matrix("`ses'",sqrt(st_matrix("`ses'")))
	mat `results'=`beta'',`ses''
	mat colnames `results'=coeff se
	ereturn post `beta' `varcov',esample(`touse')
	display _n "Abadie-kappa Weighted Least Squares" _n
	matlist `results'

	restore
end


* Mata programs
mata:
	pointer vector wls(string scalar varlist)
{ // does weighted least squares
	
	vars=tokens(varlist)
	y=st_data(.,vars[1])
	kappa=st_data(.,vars[2])
	xs=st_data(.,vars[3..length(vars)])
	xs=xs,J(rows(xs),1,1)
	xkx=cross(xs,kappa,xs)
	xky=cross(xs,kappa,y)
	xkxinv=invsym(xkx)
	beta=xkxinv*xky
	u=y-xs*beta
	u2k=u:^2:*kappa
	varcov=xkxinv*cross(xs,u2k,xs)*xkxinv
	return((&beta,&varcov))
}

end























